"""
异常值处理案例
演示各种异常值检测和处理方法
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.covariance import EllipticEnvelope
from sklearn.preprocessing import StandardScaler, RobustScaler
from sklearn.cluster import DBSCAN
import warnings
warnings.filterwarnings('ignore')

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 80)
print("异常值处理案例")
print("=" * 80)

# ============================================================================
# 步骤1：创建带异常值的示例数据
# ============================================================================

print("\n" + "=" * 80)
print("步骤1：创建带异常值的示例数据")
print("=" * 80)

np.random.seed(42)

# 创建正常数据
n_samples = 1000
n_outliers = 50

# 正态分布数据
normal_data = {
    'age': np.random.normal(40, 10, n_samples),
    'income': np.random.normal(50000, 15000, n_samples),
    'spending': np.random.normal(30000, 10000, n_samples),
    'credit_score': np.random.normal(700, 50, n_samples)
}

df = pd.DataFrame(normal_data)

# 添加异常值
# 年龄异常值
df.loc[np.random.choice(df.index, 20, replace=False), 'age'] = np.random.uniform(100, 150, 20)

# 收入异常值（极高和极低）
outlier_indices = np.random.choice(df.index, 30, replace=False)
df.loc[outlier_indices[:15], 'income'] = np.random.uniform(200000, 500000, 15)
df.loc[outlier_indices[15:], 'income'] = np.random.uniform(-10000, 0, 15)

# 消费异常值
df.loc[np.random.choice(df.index, 25, replace=False), 'spending'] = np.random.uniform(100000, 200000, 25)

# 信用分数异常值
df.loc[np.random.choice(df.index, 15, replace=False), 'credit_score'] = np.random.uniform(200, 300, 15)

print(f"数据集形状: {df.shape}")
print(f"\n数据统计：")
print(df.describe())

# ============================================================================
# 步骤2：异常值检测方法
# ============================================================================

print("\n" + "=" * 80)
print("步骤2：异常值检测方法")
print("=" * 80)

# 2.1 统计方法 - 3σ原则
print("\n2.1 统计方法 - 3σ原则（标准差法）")
print("-" * 80)
print("原理：假设数据服从正态分布，超出μ±3σ范围的视为异常值")

outliers_3sigma = {}
for col in df.columns:
    mean = df[col].mean()
    std = df[col].std()
    lower_bound = mean - 3 * std
    upper_bound = mean + 3 * std
    
    outliers = (df[col] < lower_bound) | (df[col] > upper_bound)
    outliers_3sigma[col] = outliers
    
    print(f"\n{col}:")
    print(f"  均值: {mean:.2f}, 标准差: {std:.2f}")
    print(f"  范围: [{lower_bound:.2f}, {upper_bound:.2f}]")
    print(f"  异常值数量: {outliers.sum()} ({outliers.sum()/len(df)*100:.2f}%)")

# 2.2 IQR方法（四分位距法）
print("\n2.2 IQR方法（四分位距法）")
print("-" * 80)
print("原理：基于四分位数，对非正态分布更鲁棒")

outliers_iqr = {}
for col in df.columns:
    Q1 = df[col].quantile(0.25)
    Q3 = df[col].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    
    outliers = (df[col] < lower_bound) | (df[col] > upper_bound)
    outliers_iqr[col] = outliers
    
    print(f"\n{col}:")
    print(f"  Q1: {Q1:.2f}, Q3: {Q3:.2f}, IQR: {IQR:.2f}")
    print(f"  范围: [{lower_bound:.2f}, {upper_bound:.2f}]")
    print(f"  异常值数量: {outliers.sum()} ({outliers.sum()/len(df)*100:.2f}%)")

# 2.3 Z-Score方法
print("\n2.3 Z-Score方法")
print("-" * 80)
print("原理：计算标准化分数，|Z| > 3 视为异常值")

outliers_zscore = {}
for col in df.columns:
    z_scores = np.abs(stats.zscore(df[col]))
    outliers = z_scores > 3
    outliers_zscore[col] = outliers
    
    print(f"{col}: 异常值数量 = {outliers.sum()} ({outliers.sum()/len(df)*100:.2f}%)")

# 2.4 Isolation Forest（隔离森林）
print("\n2.4 Isolation Forest（隔离森林）")
print("-" * 80)
print("原理：基于随机森林，异常值更容易被隔离")

iso_forest = IsolationForest(contamination=0.05, random_state=42)
outliers_iso = iso_forest.fit_predict(df)
outliers_iso = outliers_iso == -1  # -1表示异常值

print(f"检测到的异常值数量: {outliers_iso.sum()} ({outliers_iso.sum()/len(df)*100:.2f}%)")
print(f"异常分数范围: [{iso_forest.score_samples(df).min():.3f}, {iso_forest.score_samples(df).max():.3f}]")

# 2.5 LOF（局部异常因子）
print("\n2.5 LOF（局部异常因子）")
print("-" * 80)
print("原理：基于密度，检测局部异常值")

lof = LocalOutlierFactor(n_neighbors=20, contamination=0.05)
outliers_lof = lof.fit_predict(df)
outliers_lof = outliers_lof == -1

print(f"检测到的异常值数量: {outliers_lof.sum()} ({outliers_lof.sum()/len(df)*100:.2f}%)")

# 2.6 椭圆包络（Elliptic Envelope）
print("\n2.6 椭圆包络（Elliptic Envelope）")
print("-" * 80)
print("原理：假设数据服从多元正态分布")

elliptic = EllipticEnvelope(contamination=0.05, random_state=42)
outliers_elliptic = elliptic.fit_predict(df)
outliers_elliptic = outliers_elliptic == -1

print(f"检测到的异常值数量: {outliers_elliptic.sum()} ({outliers_elliptic.sum()/len(df)*100:.2f}%)")

# ============================================================================
# 步骤3：可视化异常值
# ============================================================================

print("\n" + "=" * 80)
print("步骤3：可视化异常值")
print("=" * 80)

fig, axes = plt.subplots(3, 3, figsize=(18, 15))
fig.suptitle('异常值检测方法对比', fontsize=16, fontweight='bold')

# 选择两个特征进行可视化
feature1, feature2 = 'income', 'spending'

# 3σ方法
outliers_combined_3sigma = outliers_3sigma[feature1] | outliers_3sigma[feature2]
axes[0, 0].scatter(df.loc[~outliers_combined_3sigma, feature1], 
                   df.loc[~outliers_combined_3sigma, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[0, 0].scatter(df.loc[outliers_combined_3sigma, feature1],
                   df.loc[outliers_combined_3sigma, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[0, 0].set_xlabel(feature1)
axes[0, 0].set_ylabel(feature2)
axes[0, 0].set_title('3σ方法', fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# IQR方法
outliers_combined_iqr = outliers_iqr[feature1] | outliers_iqr[feature2]
axes[0, 1].scatter(df.loc[~outliers_combined_iqr, feature1],
                   df.loc[~outliers_combined_iqr, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[0, 1].scatter(df.loc[outliers_combined_iqr, feature1],
                   df.loc[outliers_combined_iqr, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[0, 1].set_xlabel(feature1)
axes[0, 1].set_ylabel(feature2)
axes[0, 1].set_title('IQR方法', fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Z-Score方法
outliers_combined_zscore = outliers_zscore[feature1] | outliers_zscore[feature2]
axes[0, 2].scatter(df.loc[~outliers_combined_zscore, feature1],
                   df.loc[~outliers_combined_zscore, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[0, 2].scatter(df.loc[outliers_combined_zscore, feature1],
                   df.loc[outliers_combined_zscore, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[0, 2].set_xlabel(feature1)
axes[0, 2].set_ylabel(feature2)
axes[0, 2].set_title('Z-Score方法', fontweight='bold')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# Isolation Forest
axes[1, 0].scatter(df.loc[~outliers_iso, feature1],
                   df.loc[~outliers_iso, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[1, 0].scatter(df.loc[outliers_iso, feature1],
                   df.loc[outliers_iso, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[1, 0].set_xlabel(feature1)
axes[1, 0].set_ylabel(feature2)
axes[1, 0].set_title('Isolation Forest', fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# LOF
axes[1, 1].scatter(df.loc[~outliers_lof, feature1],
                   df.loc[~outliers_lof, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[1, 1].scatter(df.loc[outliers_lof, feature1],
                   df.loc[outliers_lof, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[1, 1].set_xlabel(feature1)
axes[1, 1].set_ylabel(feature2)
axes[1, 1].set_title('LOF（局部异常因子）', fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# Elliptic Envelope
axes[1, 2].scatter(df.loc[~outliers_elliptic, feature1],
                   df.loc[~outliers_elliptic, feature2],
                   alpha=0.6, s=20, label='正常值')
axes[1, 2].scatter(df.loc[outliers_elliptic, feature1],
                   df.loc[outliers_elliptic, feature2],
                   alpha=0.8, s=50, c='red', marker='x', label='异常值')
axes[1, 2].set_xlabel(feature1)
axes[1, 2].set_ylabel(feature2)
axes[1, 2].set_title('椭圆包络', fontweight='bold')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

# 箱线图 - income
axes[2, 0].boxplot([df['income']], labels=['income'])
axes[2, 0].set_ylabel('值')
axes[2, 0].set_title('箱线图 - Income', fontweight='bold')
axes[2, 0].grid(True, alpha=0.3, axis='y')

# 箱线图 - spending
axes[2, 1].boxplot([df['spending']], labels=['spending'])
axes[2, 1].set_ylabel('值')
axes[2, 1].set_title('箱线图 - Spending', fontweight='bold')
axes[2, 1].grid(True, alpha=0.3, axis='y')

# 检测方法对比
methods = ['3σ', 'IQR', 'Z-Score', 'Iso Forest', 'LOF', 'Elliptic']
counts = [
    outliers_combined_3sigma.sum(),
    outliers_combined_iqr.sum(),
    outliers_combined_zscore.sum(),
    outliers_iso.sum(),
    outliers_lof.sum(),
    outliers_elliptic.sum()
]
axes[2, 2].barh(methods, counts, color='coral')
axes[2, 2].set_xlabel('异常值数量')
axes[2, 2].set_title('检测方法对比', fontweight='bold')
axes[2, 2].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/0-数据清洗/异常值检测对比.png',
            dpi=300, bbox_inches='tight')
print("异常值检测对比图已保存")

# ============================================================================
# 步骤4：异常值处理方法
# ============================================================================

print("\n" + "=" * 80)
print("步骤4：异常值处理方法")
print("=" * 80)

# 4.1 删除异常值
print("\n4.1 删除异常值")
print("-" * 80)

df_removed = df[~outliers_iso].copy()
print(f"原始数据: {len(df)} 行")
print(f"删除后: {len(df_removed)} 行")
print(f"删除比例: {(len(df) - len(df_removed))/len(df)*100:.2f}%")

# 4.2 替换为边界值（Capping）
print("\n4.2 替换为边界值（Capping）")
print("-" * 80)

df_capped = df.copy()
for col in df.columns:
    Q1 = df[col].quantile(0.25)
    Q3 = df[col].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    
    # 替换为边界值
    df_capped[col] = df_capped[col].clip(lower=lower_bound, upper=upper_bound)
    
    print(f"{col}: 边界 [{lower_bound:.2f}, {upper_bound:.2f}]")

# 4.3 Winsorization（百分位数替换）
print("\n4.3 Winsorization（百分位数替换）")
print("-" * 80)

df_winsor = df.copy()
for col in df.columns:
    lower_percentile = df[col].quantile(0.05)
    upper_percentile = df[col].quantile(0.95)
    
    df_winsor[col] = df_winsor[col].clip(lower=lower_percentile, upper=upper_percentile)
    
    print(f"{col}: 5%-95%百分位 [{lower_percentile:.2f}, {upper_percentile:.2f}]")

# 4.4 对数变换
print("\n4.4 对数变换")
print("-" * 80)

df_log = df.copy()
for col in df.columns:
    # 确保所有值为正
    min_val = df[col].min()
    if min_val <= 0:
        df_log[col] = np.log1p(df[col] - min_val + 1)
    else:
        df_log[col] = np.log1p(df[col])
    
    print(f"{col}: 对数变换完成")

# 4.5 Box-Cox变换
print("\n4.5 Box-Cox变换")
print("-" * 80)

df_boxcox = df.copy()
for col in df.columns:
    # Box-Cox要求数据为正
    if (df[col] > 0).all():
        df_boxcox[col], lambda_param = stats.boxcox(df[col])
        print(f"{col}: λ = {lambda_param:.4f}")
    else:
        print(f"{col}: 跳过（包含非正值）")

# 4.6 使用鲁棒缩放器
print("\n4.6 使用鲁棒缩放器（RobustScaler）")
print("-" * 80)

robust_scaler = RobustScaler()
df_robust = pd.DataFrame(
    robust_scaler.fit_transform(df),
    columns=df.columns
)
print("RobustScaler使用中位数和IQR进行缩放，对异常值鲁棒")

# ============================================================================
# 步骤5：对比不同处理方法的效果
# ============================================================================

print("\n" + "=" * 80)
print("步骤5：对比不同处理方法的效果")
print("=" * 80)

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
fig.suptitle('异常值处理方法对比', fontsize=16, fontweight='bold')

feature = 'income'

# 原始数据
axes[0, 0].hist(df[feature], bins=50, alpha=0.7, edgecolor='black')
axes[0, 0].set_title('原始数据', fontweight='bold')
axes[0, 0].set_xlabel(feature)
axes[0, 0].set_ylabel('频数')
axes[0, 0].grid(True, alpha=0.3)

# 删除异常值
axes[0, 1].hist(df_removed[feature], bins=50, alpha=0.7, edgecolor='black', color='orange')
axes[0, 1].set_title('删除异常值', fontweight='bold')
axes[0, 1].set_xlabel(feature)
axes[0, 1].set_ylabel('频数')
axes[0, 1].grid(True, alpha=0.3)

# 边界值替换
axes[0, 2].hist(df_capped[feature], bins=50, alpha=0.7, edgecolor='black', color='green')
axes[0, 2].set_title('边界值替换', fontweight='bold')
axes[0, 2].set_xlabel(feature)
axes[0, 2].set_ylabel('频数')
axes[0, 2].grid(True, alpha=0.3)

# Winsorization
axes[1, 0].hist(df_winsor[feature], bins=50, alpha=0.7, edgecolor='black', color='purple')
axes[1, 0].set_title('Winsorization', fontweight='bold')
axes[1, 0].set_xlabel(feature)
axes[1, 0].set_ylabel('频数')
axes[1, 0].grid(True, alpha=0.3)

# 对数变换
axes[1, 1].hist(df_log[feature], bins=50, alpha=0.7, edgecolor='black', color='red')
axes[1, 1].set_title('对数变换', fontweight='bold')
axes[1, 1].set_xlabel(f'log({feature})')
axes[1, 1].set_ylabel('频数')
axes[1, 1].grid(True, alpha=0.3)

# 鲁棒缩放
axes[1, 2].hist(df_robust[feature], bins=50, alpha=0.7, edgecolor='black', color='brown')
axes[1, 2].set_title('鲁棒缩放', fontweight='bold')
axes[1, 2].set_xlabel(feature)
axes[1, 2].set_ylabel('频数')
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/Users/binming/Desktop/CodeBase/python/机器学习/特征工程/0-数据清洗/异常值处理对比.png',
            dpi=300, bbox_inches='tight')
print("异常值处理对比图已保存")

# ============================================================================
# 步骤6：统计对比
# ============================================================================

print("\n" + "=" * 80)
print("步骤6：统计对比")
print("=" * 80)

comparison = pd.DataFrame({
    '原始数据': df[feature].describe(),
    '删除异常值': df_removed[feature].describe(),
    '边界值替换': df_capped[feature].describe(),
    'Winsorization': df_winsor[feature].describe(),
    '对数变换': df_log[feature].describe(),
    '鲁棒缩放': df_robust[feature].describe()
})

print(f"\n{feature} 特征的统计对比：")
print(comparison.round(2))

# ============================================================================
# 步骤7：最佳实践建议
# ============================================================================

print("\n" + "=" * 80)
print("步骤7：异常值处理最佳实践")
print("=" * 80)

print("""
┌─────────────────┬──────────────────────┬─────────────────────┬──────────────────┐
│   检测方法      │      适用场景        │        优点         │      缺点        │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 3σ/Z-Score      │ 正态分布数据         │ 简单直观            │ 对分布敏感       │
│                 │                      │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ IQR             │ 任意分布             │ 鲁棒性好            │ 可能过于严格     │
│                 │ 偏态分布             │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ Isolation Forest│ 高维数据             │ 无需假设分布        │ 需要调参         │
│                 │ 复杂异常模式         │ 效果好              │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ LOF             │ 局部异常             │ 检测局部异常        │ 计算成本高       │
│                 │ 密度变化大           │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 椭圆包络        │ 多元正态分布         │ 考虑特征相关性      │ 假设较强         │
│                 │                      │                     │                  │
└─────────────────┴──────────────────────┴─────────────────────┴──────────────────┘

┌─────────────────┬──────────────────────┬─────────────────────┬──────────────────┐
│   处理方法      │      适用场景        │        优点         │      缺点        │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 删除            │ 异常值<5%            │ 简单直接            │ 损失数据         │
│                 │ 确认是错误数据       │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 边界值替换      │ 保留数据量           │ 保留样本数          │ 改变分布         │
│ Capping         │ 线性模型             │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ Winsorization   │ 回归分析             │ 减少极端影响        │ 信息损失         │
│                 │                      │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 对数变换        │ 偏态分布             │ 改善分布            │ 解释性降低       │
│                 │ 右偏数据             │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ Box-Cox变换     │ 使数据正态化         │ 自动选择最优变换    │ 需要正值         │
│                 │                      │                     │                  │
├─────────────────┼──────────────────────┼─────────────────────┼──────────────────┤
│ 鲁棒缩放        │ 包含异常值的数据     │ 对异常值不敏感      │ 不移除异常值     │
│ RobustScaler    │                      │                     │                  │
└─────────────────┴──────────────────────┴─────────────────────┴──────────────────┘

选择策略：

1. **根据数据来源判断**
   - 测量/录入错误 → 删除或修正
   - 真实极端值 → 保留或变换
   - 不确定 → 使用鲁棒方法

2. **根据模型类型选择**
   - 线性模型 → 必须处理（变换或删除）
   - 树模型 → 可以保留
   - 神经网络 → 建议处理

3. **根据业务场景**
   - 欺诈检测 → 异常值是关键，保留
   - 预测建模 → 需要处理
   - 探索分析 → 分别分析

4. **组合使用**
   - 先用统计方法快速筛选
   - 再用机器学习方法精细检测
   - 结合可视化人工审查

关键要点：
✓ 先理解异常值的业务含义
✓ 可视化是最重要的检测手段
✓ 不要盲目删除，可能是重要信息
✓ 记录处理决策，便于复现
✓ 对比处理前后的模型效果
✓ 在训练集和测试集上一致处理
""")

print("\n" + "=" * 80)
print("案例完成！")
print("=" * 80)
